# required packages
p_needed <- c("tidyverse")
packages <- rownames(installed.packages())
p_to_install <- p_needed[!(p_needed %in% packages)]
if (length(p_to_install) > 0) {
  install.packages(p_to_install)
}
lapply(p_needed, require, character.only = TRUE)

# set working directory
setwd("YOUR WORKING DIRECTORY")

# load data
df <- read_delim("tables and figures/fig1_data.csv")
df$factor <- factor(as.factor(df$factor),levels(as.factor(df$factor))[c(2,3,1)])

# Figure 1
fig1 <- ggplot(data = df) + 
  geom_point(aes(factor, AME), 
             size = 4) +
  geom_linerange(aes(x = factor, 
                     ymin = lower, 
                     ymax = upper), 
                 linewidth = 1.1) +
  geom_hline(yintercept = 0, 
             linetype = "dashed") +
  scale_y_continuous(name = "Marginal effects at representative values",
                     limits = c(-1,1), 
                     breaks = c(-1,-0.5,0,0.5,1)) +
  scale_x_discrete(name = "\n Model 3",
                   labels = c("Quantitative\nequitability","Quantitative\nproportionality","Qualitative\nproportionality","coalition\nagreement")) +
  theme_classic(base_size = 40) 

tiff("tables and figures/fig1_tiff.tiff", 
     height = 4000, 
     width = 6400,
     units = "px", 
     res = 300,
     compression = "lzw")
fig1
dev.off()

# load data
df <- read_delim("tables and figures/fig2_data.csv")

fig2a <- ggplot(data = df %>% filter(round(quant_equitability,1) == 0.1), aes(quant_proportionality, qual_proportionality)) + 
  # geom_density_2d(data = df, aes(quant_prop, qual_prop), 
  #                color = "gray50", bins = 20) +
  geom_contour_filled(aes(z = predprob), 
                      bins = 25, 
                      breaks = seq(0,1,0.2),
                      alpha = 0.9) +
  scale_fill_manual(values = gray.colors(5, start = 0.2, end = 0.8, gamma = 1, rev = TRUE)) +
  scale_y_continuous(name = "Qualitative proportionality") +
  scale_x_continuous(name = "Quantative proportionality") +
  labs(subtitle = "equitability = 0.10") +
  geom_point(aes(x = 0,
                 y = 0),
                 size = 12) +
  annotate(geom = "text",
           size = 16,
           x = 0.075, 
           y = 0.005, 
           label = "22%") +
  theme_classic(base_size = 40) +
  theme(legend.position = "none")

fig2b <- ggplot(data = df %>% filter(round(quant_equitability, 1) == -0.1), aes(quant_proportionality, qual_proportionality)) + 
  # geom_density_2d(data = df, aes(quant_prop, qual_prop), 
  #                color = "gray50", bins = 20) +
  geom_contour_filled(aes(z = predprob), 
                      bins = 25, 
                      breaks = seq(0,1,0.2),
                      alpha = 0.9) +
  scale_fill_manual(values = gray.colors(5, start = 0.2, 
                                         end = 0.8, 
                                         gamma = 1, rev = TRUE)) +
  scale_y_continuous(name = "Qualitative proportionality") +
  scale_x_continuous(name = "Quantative proportionality") +
  labs(subtitle = "equitability = -0.10") +
  geom_point(aes(x = 0,
                 y = 0),
             color = "white",    
             size = 12) +
  annotate(geom = "text",
           size = 16,
           color = "white",
           x = 0.075, 
           y = 0.005, 
           label = "78%") +
  geom_point(aes(x = 0.0,
                 y = 0.2),
             color = "white",
             size = 12) +
  annotate(geom = "text",
           size = 16,
           color = "white",
           x = 0.075, 
           y = 0.205, 
           label = "85%") +
  geom_point(aes(x = 0.2,
                 y = 0.0),
             color = "white",
             size = 12) +
  annotate(geom = "text",
           size = 16,
           color = "white",
           x = 0.275, 
           y = 0.005, 
           label = "92%") +
  theme_classic(base_size = 40) +
  theme(legend.position = "none")

tiff("tables and figures/fig2a_tiff.tiff", 
     height = 4000, 
     width = 6400,
     units = "px", 
     res = 300,
     compression = "lzw")
fig2a
dev.off()

tiff("tables and figures/fig2b_tiff.tiff", 
     height = 4000, 
     width = 6400,
     units = "px", 
     res = 300,
     compression = "lzw")
fig2b
dev.off()
